{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Active subspaces for sensitivity analysis\n",
    "\n",
    "A standard method to analyse dynamical systems such as models of neural dynamics is to use a sensitivity analysis. We can use the posterior obtained with `sbi`, to perform such analyses."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Main syntax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sbi.analysis import ActiveSubspace\n",
    "\n",
    "sensitivity = ActiveSubspace(posterior.set_default_x(x_o))\n",
    "e_vals, e_vecs = sensitivity.find_directions(posterior_log_prob_as_property=True)\n",
    "projected_data = sensitivity.project(theta_project, num_dimensions=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example and further explanation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "from torch.distributions import MultivariateNormal\n",
    "\n",
    "from sbi.analysis import ActiveSubspace, pairplot\n",
    "from sbi.simulators import linear_gaussian\n",
    "from sbi.inference import simulate_for_sbi, infer\n",
    "\n",
    "_ = torch.manual_seed(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's define a simple Gaussian toy example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c208f1abe76b47a29526122f41116bfb",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Running 2000 simulations.:   0%|          | 0/2000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Neural network successfully converged after 117 epochs."
     ]
    }
   ],
   "source": [
    "prior = MultivariateNormal(0.0 * torch.ones(2), 2 * torch.eye(2))\n",
    "\n",
    "\n",
    "def simulator(theta):\n",
    "    return linear_gaussian(\n",
    "        theta, -0.8 * torch.ones(2), torch.tensor([[1.0, 0.98], [0.98, 1.0]])\n",
    "    )\n",
    "\n",
    "\n",
    "posterior = infer(simulator, prior, num_simulations=2000, method=\"SNPE\").set_default_x(\n",
    "    torch.zeros(2)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "03c28be227a245129c8366307aff09b4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Drawing 2000 posterior samples:   0%|          | 0/2000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "posterior_samples = posterior.sample((2000,))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAO0AAAERCAYAAACXaiBZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAcMElEQVR4nO3da4xk513n8e//Obe6dE/3zPg2sZM4ZNloJZQ4gBLEBK0SiBSCULwCVgS0aKW8QeIFkZYXFuyiCAnhVwsvdiWEslmzEloQArSswkWRjDYkXlg73gTH9tpyQog9TjwzPX2p7qo61/++eE5117S75+Ku26n6f6TRVFV3n/OMZ35+Lue5iKpijGkON+8CGGPujoXWmIax0BrTMBZaYxrGQmtMw1hojWmY8DZft+dBi0kmebGPup+xv+cF9IXqj0/8e7aa1piGsdAa0zAWWmMa5nZ92pVz+fEnubIz4MHNNl9+7CPzLo4xb7LyNe3lx5/k4cc+z+XHnwTgys6Abz3+E1zZGcy5ZMacbOVDayE1TbPyoTWmaSy0xjSMDUQZcxckDCEIoCzRojj6gguQKHzz51NgNW3twc02Dz/2eR7cbM+7KGaRiUOCAMQd+1hO/HwarKat2eMdcye0yNGyBK1u/rws0WH6ps+nwWraU4xq3tGjIGMAUIWq9L/fyedTYDXtKUY178OPfX7OJTGLTJIElyRollENhzPp21pNa8wZSBBAuwVR5N87ORqsmpKVqmlHUxQBm6ZobmkUPM0L3+wdfT6qWYvC92PLEnb3Tq9VRZA4nmjNu1KhHc1+Amv2mtsIAh/csrxpbEnCENotJMshy9Asp8rT069TjzYrgIXWmOnRvICy9KGLYl+jVqX/LE2RbhfuvYDb26fc3jmskbX+OoBEsb9WlqHV5AaorE9rzEmqo+asxBHi/CYSWlY+oHFEudmBThsJQ9+3BVC96edw4t+PNbHPympaY+6AJAnSbiOB8xMoBkPCNyrIcqSVIK0ERKCs0KJARHwTOwioKgWtTu7TiiBhdPrXT7DUoZ3k2lhbZ7vaJI4hjg7fa5qhvf3DzyWKIAqhKCHLjn6wEKSMfCBPCaVEIVQWWuBo4GkSg06TvJZZHBKGIA4tclD1NWocI0nsv1bTCxsUG22CgxS3P0Q3Al/r5gWS5lAUaJr5GVHl0ciVFgWa1bOoTqJ6133epQ7tJNic5CVXN2H91MQSlyT+uet6F01iH8q8IL/YZXBfTPu6I6yg6kSUnYiwlxLc2Ic8R/v9N11eiwLNsxNufPP33A0L7W1YU3gFiNRN1ADptKHTRpOYKg6hHaFOKDp+oKmMHK4ToZF/L3mJ9geQ5f5aZYmO17Sn1bBnYKE1Bo6awu0WVaeFJgEaOMpuRNEKqGI/elxFjrJz1LeVvET3Dw7f+9HlW9esZ2WhNSvNra8hSeJrSFUoSlx/SHb+PNlmRLbmKNpC3KuIDirCQUnQz5G8QvISGaQQx75fWte2433hEa10Yo99LLRmdYkga1202/aPborSN3P7A/K1ezm4LyDbEIo2BJnQvlYSHuTIIMcdDND+wM90iiM/EJVn/pltkrz5XkWBphZaY27PBchogsPY+8O5xe2EqhNDHCKqVK2QKnRoKMQHStERaEMZC/ma/54wL9FOC8KjRQFOBB0M/Ej0qH87boLrbFc2tDYqvBpcHIFz9eiw+vdBgLRbSBRRdhOKbt1HFejfF5N3hWSvIt4rSTcEEMoE0g2HK0LcoIA4gG58eJ9AFent+8c71qe9e+MTIU5jo8KrQYvi5i1ggsA/hw1DCEMkL3FpSb4RU7aEMhGqEPKOowqFeE+Je0oZC0VHqPb8gJTkFVKWSFr4xQOD4cxOq1vK0I6v5gGrVVfZTc9AR8vk4giS2Ac4KwhU6b+txXDTUUWggZBFgAoXXkpp/eMWu++/n72HA5KdUWhL36/d61Fs3Ti9LzsFSxna46xWNRLF/llsHCFxTHH/JsVaTL4eUCaOg/sd+Rq4HKQCdYDAzrtjoksPULSEYACu/n+A5AUMU6gU12qho4UCM9gjaiVCa1acCG6t6wPbStAkZu/dXYYXHcN7oGgrZbdEQyW6ERAeCBpBFSj77y6JNlKCF7tsfKMiGPpQSl5Q7fV8U3vjHPQHVL3eTP44Flqz1Fyn4/uvFzchicnu6VC0Q4qWoAJVqFSRou2SoFVS7TsYCNEeBEMhSEPynYD2VaW1UxLv5QQHGaRZXbMqUpaQnzBiPCUWWrO8XIA7v4l2WmQPbpKvhex+T0jehbBfN4MD0FhpbwzZ6A747u5Fqr6wdkVZ+/aAohtStgKS7YygN8TtD+tadR9N05kNPo2z0Jrl4I6emYoTv+NEHKGdFtpJKDoBRceRdyE/pwwuVWikECo4pSwdvUGLoO8I+4Kr5w9H+znRXkZwkCH9oR8lHg6nforArVhoTfOJHD6PBfyocBIjSUJxcY2iEzE8H5CvCenFivJCwb9+/zN8oPtN/uvrl/n2zib7e22y3YTOlqN1XXG5orEjvLIDV7f8QBOgw3Tqz2Fvx0JrloNz/pFOEPjQrnXRdkJ6ISE7F7D/kCNfU/SBlM31AYMy4qXhJa73u/T7CdoPcUNHkPkR4vCgJNxNkYMBZZrWR37ITEaHb8dCa5aHc/UIcYvigU3yczE7/ywiPQ/B+3Z5eGOXRy68xkYw4H9d+16+dPA9bL++QbDvSPqCy4RoXwkyJbnWR/7xCuVg6GvWKPZ7Pi0AC61ptlHtCogI0u34tbBRQBUKVQRlorxjY5fv23yd9WBIhbCfxxwMEoJeQLwnhAcQZBAOwOUK5fFjP6p6O9V5DD3dzEJrmkvE7zRRN42JI4r7N6na4eEa2CqCsqX8q0v/l490X+Z/9N7LK/372Nrrkm23OPe60LruH+cEgxINHepAVCFJkCxD83oHijkOPo1bqtDeyZxj03DipxGi6ucUR5HfqymKkU6LYi2m6Ab07w0oWkLRUTRSelWLa2WbK+l53hicIxtEvg87UKJBRdgvCYYlaAEKMkjRBei/nmSpQnt8zrFZMiK4dtsfx5GmfoldEiOtFtXFcxRrCf0HIoYXHHs/NKDdyUi32yDwld138urwAv/n6jvZ7nUI3kiIdoX2VklrKyfsZciwQN7YotrZpYpCEFmY2nXcUoXWrICqguP9ShE0dFSxY3jekW5C0sppxzn9MEEr4bXeJlvDLtd31ih7EZ1tId6FcKhIpciwwPWHVMN64EkrJAwXog97nIXWNIeqP07yFEU7YPdflLgLKZ2wJC0CXKhUufD6N+5FMqF11REdwPmXMpIbKWXiH+W4rR3Ka9cPN2JbpD7scRZa0xyj3fjB14Rx7EeLu22G93UYXAzRTk6rlXP/+j4O5UrlGEqEXI0JD4RkR4n2IdrL/ZTEKvaL5LNsYUN6nIXWNIbEMcG99wB+h3/ptEjfdS+DeyPe+CGhXC/obAy4Z+2A//Cu/8lFN+A/X/swL2w/wPWvdei+rqy/mh7tVTxMcddyKCuqsR0VF52F1jSPiJ9XnMQUncDvMtGqcO2C9XbKuWRIr2qRa8DV4Rq7gxZBBkGuuKxEUr+Bm2YZOkzrvYonvz/xtFhoTWNoXlDt7iGtFlzYoDzf4eD+kHxdCPYdRRTyLy+9wnow5N8//yh7vTbBqy3CfaGz5weUNHJoFKD9AeXu3tG0RF28AafTWGjN4qtnPfktXfw5O1UrouhEpOeFousnUEhS0nF+Mv/udhe5EZFsCWEf4gMlHFa4tPQbjE/4+MlZstC+RXaK3uxIHBPccxHttskvnaMKHWXLsX8pJP3gPhfOHfDPO/t0woztosPrgw2Sbya0rypr3y0I+xWtKz2kP0S3d9BhSpXe4vT2BbcUoZ3lTKjxTeLsFL3ZkDCEJEbbMXknpEqEvOPI14RuO2UjGXIxOcCJ8vpgg1d7m0T7EO8r0V5J2M+RvQO036c6GMx9ad1ZLUVoZzkTymrVGXIBLo6QTofqXOdwI7Z0Xei9C4pOBXsd0tz/M+7nMa9//X7iHcfm6xXxbknynT2k16e6sY1mGQQBrtXyNW2D+rHjliK0ZgnVJ9lJu4102+TrCflaSBlD2RKKbkXVrnAKeR6wM2yzP0xIrjtaN5S4VxIOSt8kHgwOn8OOVgQ1mYXWLByJYtxaF1nrUt6/yeCBDlvfF6ICroSiDdVaiWsXnFsfkOYh1166h3jbcf9Xc+KdjHBrH0kzqmtbN02c0Czz+zo1tJYFC61ZNC443J9YWzFlNyLvOvJ1H7JgIFSR39cJoKgceRYSbztaNyDZGhLs9JG9fTTPqYbpzaPEDQ7rSKNDa0vxlovrdHAb5wDQqgIRyiRAKqV1ze/+X0XgCsHthrg8gutt1gbQ2qqIexXhd7bR/QPKfh/KsrGPdW6l0aG1pXhLJgj8cR2qSJajgUOdIArhQCnUH4QF/ujJoC903lDCoRLvlYT7ue+/pqkfdFqCWvUkjQ6taah6VPhwJc1olFgE0gxaCdW9m5RrCWUiSAnt6xXDC47e91ZIJnRfc0T7SnurINoviF95w29tWh8z6ZLE76C4hOG10JqZEyf1FjF+y1MJAogiiEK/4iZwaBJRxQEaCCiEwwopHRoqUvgN2OKe+h0n9jPK637A6XD7GeeQqprLZuLTZqE1M6dliQ5HO08kuPU15Nw65fkuw/s6hMOSaHuIaEjWdSAglQ/vxgsh4YFy4f8NCA4y3HYPHaaUuT+io0rTm4+2XLJaFiy0Zh5UQUtwoZ9PHNUjxZ2YdDNAe0J0o57g7+pfkSCVbybHvYrwmp+WWO3u+XN0RgNOo2svMQutmRsJQ3+KXZ7DtW1C54g3I6RS8gttpFDOv3RAthGz984IV0D7ekG0l8Fuzw84DQYLuSXMNLnbf4sxUxIESOR3otB+Hxmmfs9hhaLlZy4F13tEe35wSSoIezlBL63P1EkbvVrnrbKa1syc63SQ9TXYWCe/uEa4O0CubgH4xzaRq/ctdpQX1tBAWH+1IByUhNd9s7gcDBu1cH2SrKY9o9Gqn8uPPznvojSGxBHSaVOtt8g2Y6pu4pvKVb2zRFH3ZwOh7PiaONkaEt3oI/t9tD9Ai3zlatgRq2nPaLTqx5boHZEwhCBA82NN19EURYAsx+2nJIHDHaSoKkQhZTdCnRBkFa6oCA5yJC+RgyFSlEs5Gny3LLRm8sT54JblTYfMjXafQJxfcZNmBAcB1BMiCNzhtEWXVUhW+Q3E08zXsBZYoKGhtTnHDeMCJAhw3TbSbvvFAFGIjpbJOYfEMQwz2q/uoSLgQAYZbO9CWfkjOio/w0lVV7rGbWRobc5xs0gQ+GMikwRaCRqFfvbT4TcIhAHkBRz0cWG9yqc/oNy6AeJwrQSq6pabla+KRobWLDYtS0gr30yOYv9ZmvrphZWCE6Q+S1aTGMkLH9jAId2O3950e8fPI1YFKl/Drtjz2NNYaM3kVb4vK/VBzJrl/nnqMPVn8QAKSKfjw1uUaJ4jQYK2EyQv6u1N65CqNmb3/1mw0E7I6NHP6PVK7iVV91218DXkaPcJ8tzPenIBEji0rKAofKDjCAHIQPtDqA/B8kdZ1seAaGWhHWOhnZDxkK7q4x+J/FziUciklVBtrCFZ7pvAo+/LC9SJ79eGwdEkiTSl3Nsbu6BDotB/3UJ7yEJrJkbz4qYjNnQwwO0EaJ77/mng/NTFVgLdc37t7PXtoyM6jgdTrS97kkaEdvSIZ8Qe9Syo6ubnstXBwNeqqZ8jPOrjSrtFudEmvJpRvnH19OtZX/ZEjQht0x7xjG9ovhJ9WxcgTnwNq+rDOWoqB4F/DluWSLvln9MWpd/Lqd9Hotj/3IpOSXwrGhHaplm1qY2HAR3VqK0E6Ry1hiT284el00Y7LWRrh+K7b9Rn8yR+EMpCe8cstObM/BxjPZq+CH5qYt2HZfTIp1K/BUw9QULL+sjJFV2t81ZZaM3ZVSValUiSIKGfuaRZ5qcmBgGaplT1elnp9w83X0O18efqzIOF1pxd3af1c4Pzo89V/aOaqjqsgdVq1jOz0JozG/VpR6PEI6NaVKIYSRIfWKtZz2yhQ2ureRbceA1bh1XC8GgUOUl8E7n0hzhbDTsZCx3apj3qWTUShUgQ+PNy8lGf9mgU2a11kXPr6PYO1U5/3sVdGgsdWtMAIv60AK37rEVxNIMp94sENLcJEpNkoTVnJt0O4gK017tpvWs1GCJlefPglDkzC61568rSH7uROUQKv3pnXL1wwPqyk7UwoR2fX7wy0/8aTosCCj+3+FZfN5O1MKEdH3Ralel/jSPiz8mxKYdzZfsemzsmcez3anLBvIuy0iy0U7R0G5lXWj+DrW7/vWZqFqZ5PG58aVuTLdtqH5vNtBjmHtqTZj3ZINSCGfVlwdeyK7zn8CKYe2ht1tPikzj2p7XXqmFqg1FzNLfQ2rziBSLifz+tBq3UnzM5Yn3auZpbaFephl307Wdcu+1370/TE4OreYbapKaFMffm8SpY+AGpqrLDrRpk6qEdbwZ/+bGPWLN4Adn5OM0y9ee0J810+tbjP7GQzcRpW7rntmYu5DbNorfcZjpew5ojE/hvI5Msz0fdz1jbeAF9ofrjE/+ebxnay48/qeNN2pMcb/Ye/9yc7gzhtdCugLcUWmPM4rG5x8Y0jIXWmIax0BrTMBZaYxrmlpMrROSvgHvOeI97gOtnvMYkLFM5rqvqxyZRGNM8Ux89FpFnVPUHp3oTK4dZIdY8NqZhLLTGNMwsQvt7M7jHnbBymKVgM6KMaRhrHhvTMBZaYxpmqqEVkZ8XkX8QkedE5CkRed8073dKGT4mIi+JyCsi8tis71+X4e0i8jci8oKIPC8ivzyPcpjlMNU+rYj8MPCiqm6LyI8Dn1HVD07thm++fwC8DHwUeA14Gvikqr4wqzLU5bgEXFLVZ0VkHfgK8Oisy2GWw1RrWlV9SlW367d/Bzw0zfud4APAK6r6TVXNgD8EPjHjMqCq31HVZ+vXPeBF4MFZl8Msh1n2aT8F/OUM7wc+GK+OvX+NOYdFRB4G3g/8/TzLYZprJrsxisiH8aH90Czut6hEZA34E+DTqro37/KYZpp4TSsivyQiX61/vU1E3gt8FviEqm5N+n63cQV4+9j7h+rPZk5EInxg/0BV/3QeZTDLYdoDUe8AngR+QVWfmtqNTr9/iB+I+lF8WJ8Gfk5Vn59xOQT4feCGqn56lvc2y2faof0s8FPAP9UfFbNe4SIiHwd+BwiAz6nqb87y/nUZPgT8LfAcMDpT41dV9S9mXZZT2LS4xXPq5n02jdGAhXYRnRpamxFlTMNYaI1pGAutMQ1jp+aZpbIKx9FYaM1SGR34dvxY0fFja5oe6MY3j0XkMyLyK/Xr3xCRHzvDtT4nIldF5OuTK6GZh+MnFI7CPH6CY1NPL1yqmlZVf/2Ml3gC+E/Afzt7acw83eog74U/5Ps2GlnTisivicjLIvIl4D1jnz8hIj9dv/6WiPxWPZ3yGRH5fhH5axH5hoj84knXVdUvAjdm86cw5q1pXE0rIj8A/CzwCL78z+LXp57k26r6iIj8Nr4WvQy0gK8Dvzv1whozBY0LLfAjwJ+pah9ARP78Ft87+tpzwFq9lrUnIqmIbKrqznSLaszkNbJ5fBfS+vdq7PXofRP/h2VMI0P7ReBREWnXW7f85LwLZMwsNS609bYtfwR8Db8TxtOTuraI/HfgfwPvEZHXRORTk7q2MZPSyCZivbzuTUvsVPXfjr1+eOz1E/iBqDd97djPf3JihTRmShpX0xqz6iy0xjSMhdaYhmlkn9aY48ZX94wbzUE+/nmT2XYzBpZgu5mHH/v84WKAO7Xgy/hO3W7Galqzspq6cMD6tMY0jIXWmIax0BrTMBZaYxrGQmtMw9josTG1BX8EdMhqWmNqo83fRrs2LioLrTENY81j01jH9zJeFRZa01ij5uyqseaxMQ1joTWmYSy0xjSMhdasvKatubWBKLPyFnkixUmspjWmYSy0xjSMhdaYhrE+rWmc0zZxWxUWWtM4qzoTasSax8Y0jIXWmIax5rFpjFXvy45YaE1jrHpfdsSax8Y0jNW0xpzi+CL7RZnuaKE15hTjzfFFOjrEmsfGNIyF1piGsdAa0zAWWmMaxkJrTMPY6LFZeDYT6mYWWrPwZj0TatH3jLLQGnPMokyiOI31aY1pGAutMQ1joTXmDoz6uZcff3LeRbE+rTF3YtTPXYQ5yFbTGtMwFlpjGsZCa0zDWGiNaRgLrTENY6PHZmHZnOOTWWjNwrLdF09mzWNjGsZCa0zDWGiNaRgLrTENY6E1pmEstMbchUVY7WOPfIy5C4uw2sdqWmMaxkJrTMNY89gsHJu+eGsWWrNwbPrirVnz2JiGsZrWzNV4U3jR9xseN76h+azLbTWtmatRU3h04npTfPmxj8yt3BZaYxrGQmtMw1hojWkYG4gyC2E0sDN6bU5noTULoUkjx/NmzWNjGsZCa0zDWPPYzNRoMsVI0/uv85hkYaE1M7Vs84rnsb7WQmumZrxWbdo0xUVmoTVTM16rLsK5rsvCBqKMaRgLrTENY6E1pmGsT2vMBByfhjnNQTcLrZmJ8eeZy2g8pNMedLPQmjO7k90nVulxz3itO3o/yT+/hdac2ejRzuXHn3zTP9ZVdDygk655LbTmrp1Ws65SbXo3Jj3V0UJr7qgpd3x203jNuqo16p06barjW50xJqo62RIaY6bKntMa0zAWWmMaxkJrTMNYaI1pGBs9NojIXwH3nPEy9wDXJ1Ccs1qEckyiDNdV9WMnfcFGj81EiMgzqvqDVo7pl8Gax8Y0jIXWmIax0JpJ+b15F6C2COWYahmsT2tMw1hNa0zDWGjNmYjIz4vIP4jIcyLylIi8b07l+JiIvCQir4jIY3Mqw9tF5G9E5AUReV5Efnkq97HmsTkLEflh4EVV3RaRHwc+o6ofnHEZAuBl4KPAa8DTwCdV9YUZl+MScElVnxWRdeArwKOTLofVtOZMVPUpVd2u3/4d8NAcivEB4BVV/aaqZsAfAp+YdSFU9Tuq+mz9uge8CDw46ftYaM0kfQr4yznc90Hg1bH3rzGFsNwNEXkYeD/w95O+tk1jNBMhIh/Gh/ZD8y7LvInIGvAnwKdVdW/S17ea1tw1EfklEflq/ettIvJe4LPAJ1R1aw5FugK8fez9Q/VnMyciET6wf6CqfzqVe9hAlDkLEXkH8CTwC6r61JzKEOIHon4UH9angZ9T1ednXA4Bfh+4oaqfntp9LLTmLETks8BPAf9Uf1TMY8K+iHwc+B0gAD6nqr85hzJ8CPhb4Dmgqj/+VVX9i4nex0JrTLNYn9aYhrHQGtMwFlpjGsZCa0zDWGiNaRgLrVl4IvIZEfmV+vVviMiPvcXrzGQVzrTZNEbTKKr662f48QL4d+OrcETkC7NeDXRWVtOahSQivyYiL4vIl4D3jH3+hIj8dP36WyLyW/V0ymdE5PtF5K9F5Bsi8ovHrzmrVTjTZjWtWTgi8gPAzwKP4P+NPotfm3qSb6vqIyLy28ATwGWgBXwd+N1b3ONhprQKZ9ostGYR/QjwZ6raBxCRP7/F946+9hywVtegPRFJRWRTVXeO/8C0V+FMmzWPTdOl9e/V2OvR+zdVSrNYhTNtFlqziL4IPCoi7XrA6CcncdF6Fc5/wW+P8x8ncc15sNCahVMPFv0R8DX8ThhPT+jSl4F/A3xkbD3wxyd07ZmxVT7GNIzVtMY0jIXWmIax0BrTMBZaYxrGQmtMw1hojWkYC60xDWOhNaZh/j84kOvfMCHqPgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "_ = pairplot(posterior_samples, limits=[[-3, 3], [-3, 3]], figsize=(4, 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When performing a sensitivity analysis on this model, we would expect that there is one direction that is less sensitive (from bottom left to top right, along the vector [1, 1]) and one direction that is more sensitive (from top left to bottom right, along [1, -1]). We can recover these directions with the `ActiveSubspace` module in `sbi`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "72af7fbea54d4facb49911908cf463f7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Drawing 1000 posterior samples:   0%|          | 0/1000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sensitivity = ActiveSubspace(posterior)\n",
    "e_vals, e_vecs = sensitivity.find_directions(posterior_log_prob_as_property=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The method `.find_active()` returns eigenvalues and the corresponding eigenvectors. It does so by computing the matrix:\n",
    "\n",
    "$M = \\mathbb{E}_{p(\\theta|x_o)}[\\nabla_{\\theta}p(\\theta|x_o)^T \\nabla_{\\theta}p(\\theta|x_o)$]  \n",
    "\n",
    "It then does an eigendecomposition:\n",
    "$M = Q \\Lambda Q^{-1}$  \n",
    "\n",
    "A strong eigenvalue indicates that the gradient of the posterior density is large, i.e. the system output is sensitive to changes along the direction of the corresponding eigenvector (or `active`). The eigenvalue corresponding to the vector `[0.68, -0.73]` is much larger than the eigenvalue of `[0.73, 0.67]`. This matches the intuition we developed above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenvalues: \n",
      " tensor([2.3552e-06, 7.0754e-05]) \n",
      "\n",
      "Eigenvectors: \n",
      " tensor([[-0.7066, -0.7076],\n",
      "        [-0.7076,  0.7066]])\n"
     ]
    }
   ],
   "source": [
    "print(\"Eigenvalues: \\n\", e_vals, \"\\n\")\n",
    "print(\"Eigenvectors: \\n\", e_vecs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lastly, we can project the data into the active dimensions. In this case, we will just use one active dimension:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "projected_data = sensitivity.project(posterior_samples, num_dimensions=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some technical details\n",
    "\n",
    "- The gradients and thus the eigenvectors are computed in z-scored space. The mean and standard deviation are computed w.r.t. the prior distribution. Thus, the gradients (and thus the eigenvales) reflect changes on the scale of the prior.\n",
    "- The expected value used to compute the matrix $M$ is estimated using `1000` posterior samples. This value can be set with the `.find_active(num_monte_carlo_samples=...)` variable.\n",
    "- How does this relate to Principal Component Analysis (PCA)? In the example above, the results of PCA would be very similar. However, there are two main differences to PCA: First, PCA ignores local changes in the posterior, whereas the active subspace can change a lot (since it computes the gradient, which is a local quantity). Second, active subspaces can be used characterize the sensitivity of any other quantity w.r.t. circuit parameters. This is outlined below:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computing the sensitivity of a specific summary statistic\n",
    "\n",
    "Above, we have shown how to identify directions along which the posterior probability changes rapidly. Notably, the posterior probability reflects how consistent a specific parameter set is with **all** summary statistics, i.e. the entire $x_o$. Sometimes, we might be interested in investigating how a specific features is influenced by the parameters. This feature could be one of the values of $x_o$, but it could also be a different property.\n",
    "\n",
    "As a neuroscience example, in Deistler et al. 2021, we obtained the posterior distribution given burst durations and delays between bursts. After having obtained the posterior, we then wanted to analyse the sensitivity of metabolic cost w.r.t. circuit parameters. The framework we presented above can easily be extended to study such questions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d02a088dbde245049d060f27e51a4b65",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Running 2000 simulations.:   0%|          | 0/2000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Neural network successfully converged after 139 epochs."
     ]
    }
   ],
   "source": [
    "prior = MultivariateNormal(0.0 * torch.ones(2), 2 * torch.eye(2))\n",
    "\n",
    "\n",
    "def simulator(theta):\n",
    "    return linear_gaussian(theta, -0.8 * torch.ones(2), torch.eye(2))\n",
    "\n",
    "\n",
    "posterior = infer(simulator, prior, num_simulations=2000, method=\"SNPE\").set_default_x(\n",
    "    torch.zeros(2)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "18a9275d763e4f07907a0ed7e93591c4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Drawing 10000 posterior samples:   0%|          | 0/10000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAO0AAAERCAYAAACXaiBZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA0rklEQVR4nO29WaysWXbX+Vt7f0NEnOGOVZlZWVkkBhrTlvFAtZHIQgIPknEJ2RIgYVCjlnhB4gFLzUMKJNpCQs4n4JFGhVWohRraAiQkY5BRumVwCXfZiaHssqpc5ZpyqqzMe++ZIuIb9l79sPYXEefcc07e4Uxx7v5JR/fEcL7vi7jxj7XX2msQVSWTyawP7rIvIJPJPB5ZtJnMmpFFm8msGVm0mcyakUWbyawZWbSZzJpRfMjjeT/oaiJnebAfc39pPf6fRZCiBCeI93ZfjKgq2vX2lHS/9h2Iw29vgjji3h7a90hRgDh7XPWh2498KWWFjGq0bdGmOfOXCvDL8ReO/X/+MNFmMpfPIFaNaN8hVYVsbSIiJsRBbCGg8wZVtecDcf8AjYp4j5RVEme/OLT2/XFn/FC0a+1ahuOGADE89Ut9FPLyOHO1EHn4tjhw6d/hbueWj4mA87CwvmrPBxOTxvO5VlU0XvwiJVvazJVBigKpaxND1yPeIeOxWdgQERHwG6CRuLtnf6RqlrgsICqx7UDjchmdLC5OTMxnTQyoxpOX1mmVcJaWOIv2CK+89jpvPZjx4s0xv/bqD1/25Tx7OAcxmsi8B+9ABRmsp3cQsKVwEoqIoMPfHbWqg1gHwYqziMBZLmUfwRcWJ2dm8PPy+AhvPZjx9dc+zVsPZpd9Kc8OIkhZARAPpsvATgjodIYeTAn7B+afTmdo20JZWhAJ0BCJB1PibG4+ZlEilT2uXY/2nf2EgL9zm+KF53AbGxZMKqsUjLLltNS1PZaOjTOf9eiyXcpqcYxjX4tLS3XVdP4n852PI4s2c7WIwZaSYJY02tIYjeY/RjNXIhZBFpGF/ypObIldFhZFdunjrbqwhlJ4GB4fWPWVRZZ/NxzTu0PPAeyczi3OfRRZvf+Mi3Ke6eXxsBQG8nL4PHE+LQ/1+GVpskZHP9ziPVLXy/2tokA2Jybk6dx817JAnEcmIxORiD2+tw9tZ75uSEcQR7z/AGDh+66eCyyarCuPSVUh4xGiB2gTFq+FEIgH04eX4ye8lrPkmRbtsBQGePnVX7zkq7nmeA/xlCXi6od8sGoiUCw/olKVaFWatZ03EB14s646qlOEWaBffjGICLoSVY5Nc0rQaGW/VtKxvDeher9MWvDe9oVP8ovPudz1mRZt5oKIAW2PibCmJS7eREGM5vsdXYr2PXE6NcubrKA2rS2Fb92wQ+3u2/FHNYRoljQEZFQj4tDZDHUsoskyqhHvTHyqliSR9l0Rh9veREYjdG+f+GAH8Q43mZiP3LYgzr4sTlg9mM87Ts/vltHj9JpPXHU8As+UaE+LDL94c8zLr/4iL94cH7o9/J6Xzk/JcdZH3EKwUhRpT/WY52k0oaTnARACFAVaePubbiXDafUYKeMJ5+xvnJg/XJVmQbFA1qFrdIJUFdQVHAjateBqXApsmSXm9NWDs1WCSDqvxuVKWhzi4hNHk58p0Q7L4Vdee/2QQIGHRLl6Oy+dnwyLzhbJiq18uAcLmyyqWaGVT3A0v1CGiG1Z4e/cNkGOR6h3i31ZSemLevumPbc34bs7t0xcTWt7vEVhoq8r+zcFm+K0R5vGosaTifm0IaCzmR0rKm4ysXO07UM+rAW8HvbXtW2JO3vmc4/HMLwHg8/7FDxToh3IVvOCSNFcDQEOidYts5cSh9oepWQFVUXAordFBYVHxzVaOOKoRLqIf2AC0FFpPuisNbGWZnGlDxCiCRWQslz6yaq2rdT3JsyU2EEIltwRh4hzYWI+LgPKCeARwmE9qyYLnSLaIRx67Gl4JkWbuRi0bYnpwyorASWN5kO68RjZ2DCfd94s/Fvte5jNLHI7GUEf0P19E9W4RvqIfzAFEeLWGCK4vdnytgjSpaV2VZpVHqx2jMdmRmkISCookKKwpa136HyOtt0iW0ubBm3bheVG1b6UTvJxQ0Cn08XS/Wl82YEs2sz5oWoCdN6EMxCTX5iWq9p1iO9NBEVhfqb3ZmG9hxBN1EOAKkSkadGyQOuxCfF+ByJouYk6wfWW9qjJ4pLOL7PmYdGILCwsVbn0nVNus4ZgPnCymBrCwhcf9pRPsriHVhn+GIv8BGTRZi4NnTeo37dAUAj20/f2L9h+6f0d8A536ybUFbH0Js66QL0njAqzvMmahnGJFoJretQ5+tsT8ELxYI50AR3bElj2ZxAC7s5tnAjapWhzXUPhzRfueygK3Pb28qKToGVU21J7PrfnhWVSCHA4Ml4Utq97ZG/4ScmizVwaGgJxNj/1cZ1OceMRbG+Z1XQO9YKWHnVCrJyl9TmHekes7PFhOdxtlagX/EGH9BEtzFoLoDHC5sQypPamFmgqvJ1n3kDbwXhkkeS+R/uwrOV1HsrSotaLCz6y1zxExr1P2V25YOCROJr1lLlEQjCfLlkbV5XI5sZyaRrNPxTvbZmaGCwVTqCPiA5RZQ8KsXT0H91GglLuWN5yrDzqHb6JaCF02zWyWVF9ew+ZpwDRqDYx9QGtzffVISWyLJd+cG97s9p2qShBLcNreB3jkUWou862llb3nFf+fcinHfKTH1PM1160q1lPT8rqHm6OPD85DwVhytK2Q2JMggiWsliVyGh0OEk/KhKiPRcHDiSYZVMH3WaJ65XqPUuy6G9N0NLhQkSj0G3aR716y4oQZDJGffKRY4SyQIdsqhhtiSxigo0RhqSKgRaLLpcFlBXiwzKDqiigaxcR52G5f8inlWX3Dc2iPXsGoeb92idEY0pKOOzPadui+weLelnA/EpnYmFI+vduuTQuHFoVtDctbdF1EQmgTgi1Y/d/voUKuN5E5ToTdizTcvm5bfzNDQZj7XdmSNuhlaCFsy+NIMSNEVp6pO2RPiJ7U9xsZiLu+mVRQdejutxTHiLfqy1wjr4XgEWby8K+hB5z2zaLNnP+qII+bE20s/TEh+gEymgCGooAysJqZp0jVp72RoFEqB9EJKr5t6Vj70WPehjdU3yr+EaRCDF90ttbFRIU1yvSK35vbpZ1yFGOlq8cNmv6scc3Ja4L+LRnrPNmsX8Lye9OXzBSFNB1xPnJfvohvH+4U8cjkEWbuTSGUroFw56nc1ZCB8vtkj5A5QiTEi0d5UHKZR5q26Pi54GtN4VQCfPbjnZbKKaKBPCt4npwbcB1SnuzIBZCrG7i5wE/7ZAuoiNbJsfKrDrzgPTRrPyoWmQ2yXiEjEZIjLZ8Hpb+qsjKFs+CQdyL12rbWMBj799m0WYuD3GH6lqHLRMZfMO0zysxLtIcw6gAgWJ2eOkpUfG9Mnl7ThgX7H+8ph9DqATfguym57QR10b6UUm3IfSjAtd7Nt6GYqexAFbpiKVDh7qFoKiXVAboiCEgzrKzFufverPCIaQ9Zn84oaTvj4g2ZUyldjTiH923zaLNXD4pWiwhWBWP91bFU5WwtYF6j44r4qhAC1kshVHzWbUQ2m0TfzGNxCpFgBX6DfvpNgTfCr4rKQ4irodiatZXFOZ3SuRWiW8j0ivFNFB2wYJdhQlYvcfdvoHbnCyL6kO0VYB3sLVhOc2z+eHo8Wow6ijD6uIxyKLNXD7e2V5ojOZXhgBNA1WJbo7R0pIotBCiF3AWWJIgiEZC5ZjdNmHVO7K0kEC/ocTCBOZ6odpx1E5wveJWAkDzW45+JIweCOU0Ut1v8Pf20c0xcVQslrr9pETdBL/f4vamtl3UtsjGhLg1waeA2kKoR5fFx5Gjx5m1QWPK+RXbThn8utVAjzer6ec9sXRI5dDUy0I9hNolAQNi4oslzO8qWijVAwdRaO5EYqnEUggldBMTuetAVJnfEfoNKKdCtQexLpCbG/RbFf3IUz1oKXbnltRRelzbm5UltXMNAbc/s24ZRZH88lRoMHS4OKPi+CzazKWhUS3QpHooIrtAxMrwouJmHVQFblQQSwBbJofShO16a2TR3IUwUcLzDahQfrOmmCntLYhjJVSOUEO7LcQSiilIgPkdpb8RGL/rkaD0E0+oHd12QTcWql1Bdg+gLq28r2mXHTKKwn6fN+aDD8v98RjmDdK2j2ZxH5Es2szFkxLxF83PYlwUwA8ZUVIU0Pf4g3aRthjHBV3yXV2nxAJi7YiF0I/TkjmC6wRtPIj5s4Mllga6LQgjodtUtLD9Xd8IvlXY84SRMP1owehBoNgP9LXQbVpyRrE5RuuSWHncsFXjLe1RmtaKEZxDnbfot3dIVVqd7mpNsfOnW17nj78/kUWbuXiGxIKEpRL2S2uUyuO075GdfRjXhFsbhHFBs+1xvTL6oMepCTJU0G8I6sEFLFupcWihtNuKCynZIgjtzYh6iJMAXpFYogeCnwt+Dv0E+rHgW2eiHQvttgW6qhtjYu0JpaN0go+ROKkI4xI/9Ti1rC1JgSktrFukBaccHBwsG6ufYnmHwvqTuLaiXW0tk1kTvMcVhSXip1RG3d4AEfx+kzKgLCe53/AWcFJb3roWYgXNbSV6cHOHeqW/24GAv18gQdDnGsq6p7k/QmaeWEA/Uco9wTcsgljtphCLCnVQ7lmCRhgVxMp86F7LRRqkxNQIriqteD51z1hElrvOCvEHCzrkYJ9gcTWEZX71MVxb0Z5FzvFRcg7y+SJVZQkL3kFZopMR3e0Jftbh3/4Ar0rRjAmV0GyZYFxQXC8Uc6XzQnu3h0Kp3ypREW599wNu1HO+3H4MN3f8sZfe5RMb9/mlB99LsecIYyWMlfqep9pR1Jtw25vCdAyTd5XxPcWt+LmxFPqRwwVPsR8o962WN9YlbhDgEFTrOiukD/GwhR2GiuFtv3aVeHrN7bUV7XmQc5DPlkVj8PSvDN0iVGE6tVk+UYlVAc/fQUtHcRBwrUPFor+xsC2ewUK6mUdLJYwglop3ERFFJhZ93i7nbPgGN+kJE496E1m3bUttW0ZDP4IwsmizOrVigpXJr+Y726iSWHtcE3Cz1NoGFpVBiEBdW37zbLZsV6O6bMD+mGTRZi6PJNRFx4pkYZnNCPfu4wHX3aDfLJm+OMLPlfG3Z5ZYITX9SOhHghbY/q1CueMItdLf6pFRwIuiKmzdmKHAC6Md7pb73Nie8iAKOi2QXpg/14NA+cDjZxaoiqNIqD3Rgz+yWnW94pqIeqHbKKg666Yh/VDA3xI+uIfb2sJ95A64OfrgwaGl8JNGk7NoM2fPEFk9tm1qemzw62K0Ej1AXWn5yGWJ29xEyhKZdbjC4foCCWpVN1IkC7vc6um3rShA0iknd6dsT+aUPtCvZFt8u9mmU8/Hb+xwYzxnZzai7QtCcMQohIMJfipoqWgdianYvq8hlOYDx0IAWwVIr/h5SqjYGEHTpSiyP+zD9ivtYY8WCRy9/0P2c7NoM2fOokHb0S78qVeUDcVKrU87FuM7nNtAyw3rJLExhhBxeweIKv5GhW8CMu8Q5wiVRYuLudI5oblt/mi5Zx/8H/kDX+a7xt/h9fe/mwfzMQrE6Pja7m3erzb49HO/zXfV3+at7jY7Ycw3Znd5v93g1+/9Yap7jlhFis2OUJWoh3Ys9kXhLakjFoA66p1AeX9OHBe0d8YUBwW+65HCL0vvpnMY/NaVafbAMT4uD/u4R8iizZw5GhXhmGHOGpf1su6U5uTOWduXPiBTy+OVqMTS0T2/RSwcvlMCtjzux0IsLAMqFgIC3zy4DcCkaHEjJaqgwB/cvsft6oC5FnyrvcNcSzzKdmHdTYobLU3jYGS+cD9R5rdlEaCStH0EHLaYaj/qBB3VlhAymdherXdov5JbqdHGdab3yn55dB83izZz9pwU/RyqdsoKNx4Rm2Y51nKVwhM2a9y0s3EfXY/0Srfl2Xm5xHWw8W6PROHghYJuA+I4gkCs7UP/xbef51ubN3jlha9xa3PK26ObAPzUnd/kI36PX7j/Q/y32Sd4aXKfG37G3XKfu+U+P/CJTd6+fYMH0zFtU9Df6dmfOMtzjuYzV7sPX7KoWpcM7+hvjfGjAudkue0z+K9Dh8rj3p5H9HGzaDMXw9CmFBYDo61IQG2W7NDdsA8o4A/E6lrv3rIa1yE63C0tnUoqDkj7teqUWNoyWVSYtyVfuP8xxkWHE6VyPb/XPM+DYoMH3YR5KLhb7vFSeY8vzV/gg26Dke94fmOXvXlN6B10gmvNeqtAGCutCK4T2LO0yP5mvbgeFyOu6c337vq09ZNylI/WzZ7kw35IYXwWbeb8WfFlialxW8p6ktHIGqxtjqHpLLk+5STL83fZ+d47oEo5tbYyowcRUbVKn4KUighEAa+Ebct08kAzr/jmN1+wYoKPzZlMGva6EeOiY9aXeIn8L+Ov8cl6yq/sfDdvvPcSP/jRb/G9W2/xtQe3idOCctdT7Av9hhJGSncr0I8CEmom34F22zG/XVPMlXI/4JuAfzBFZg26t7fsGQWLsSNxOl36sBofsrCHfN5juHaivYhMqDyc6/GQNPSKYVykiPl1Q/NxVUiDrRbtSpM1cr0iwSK01mnCE50FhfqRI1Qm2vq+BafCS3OcV/p5Ab0gaQ/XAarCuOjYLucUYsf/avtRAt9hrxvhnVnEJpaUPkIViR6zsg60UKRxyNzheujrVLcbl0Or1Yk1hUsN3wY/HrAsLydIl6bHH+fDpiHZcoq1vXaiPY9MqKPk4VyPSVoKH/oYtp19sJ1H53M4mMLGhHhn25aW+5aIMPqgXfRyCps1za2CUAndBPqJ0NyO+Lnw0Tci3cTRfs+MrVHDO7/1PL6B5uMdxainKANOlD+2/S5/sP4O73Q32enH/D9vf5LdtmZSdtwczZiFknebG4zLjq1bU/bmntB4wjii48D4GxWb31JCbQGqalcZ7eiiM2QsnfnjpQ2flnlrI03KCtnasBK+soKuJezsHhuIk6p6tkSbuWAORVA/PPI5DNoaPpQWWbUlooTUjymkfNwYkTZaMCdlTklMXRaD4DooDqxKB+z+nQcTpqMK9UoYCdVGy2TUUvhIVfSUEgg4pqHioK856CqmTUUfPIUPbJYNtQsULlL6AGUkVj6ZarFc5REgyb8OWORYIFZWmL8Ybi1mdWU0srGcLtX9auopJQ6IK61UrZ8yMVoTuxPIos08MTLMfk0spqifwNC5cJEF5ZK/NzQDnxfI/tRm93TWlrTYnROrgrBREWpv7WCSYMt9pdq19MJ2w9Iab/x/I/oJ7P9PHeVWw/e88A53auv46CTiJfJ+t8XXDu5wv5lw0FR0vWdvZ2zbMh8Dt6U4lI2qY2+zpVOBIEjr6LYje2Nh/K6w+XbE9fZ6+zp1vgB0d/lFpnWFbk/MErdWOKDTmRUFlIWVIY5H9uQu9UnuulMjyVm0mSdGoyInlX6uRoujLi3syqQ5icsG39btITULF0kjKe3gouYzuhDxM0HKZXKF9W6ycjr1VukTKiAIofO8P9ukjYX5sS7w9vwmAG0scKK0vadtShNsFJwohUR6dcy6EucUqQN6UODmliq5KMBJ+7KgFgwTUo5yegtCRMuCOKmQJuCaZdKEiCzHeya/dlFTnDOiMudGDMd3EHQrE+VimoLetri6tvK12Zw4neKcQyY2gkNu3lg0Q5NRTbyxsVhiqthALZkq5Xsd/c0JB89PrJh9w4Q6fy6gdaTcbIkqFG+OcR943vrOc7xZKZsf32Wjbtk5GBOj8PLde9yqp3xl96P4eyVsWOBpq2p4YbTDl+59hPv3NilHPZvbM/bv36B+31mec2FJFjEV0auzLaHD740FoLT0zO+OKPd76vt7llyyMrNXQ1g2N08TDJ656HHm8rGMJzmU9SOyHBKNW2md2veLDCiiQ1RthGW50shbrPuiRDXXsnArGUgpU6kX1AtRbUNVvRLjsA8KXVcwBZp5iQZh2lUULkJviRN2HrNwAUcXPNp6tFp+KYkCPfhgv1sFkFl66ZdTDUR1OSIE8POAa4Itf/uHv+Q0xMN1tRotofoEsmgz589gkZ1YuV1RmGXte+LOLjIZI/WWJVHUJtgwXg7gCiNPc7tEZVkQUB0ofQ/dphWvb/2+Tc2bfcQTayXc7nCFZUk5ga4taJuCuFsiQXi32uZePUHaw+KY9yU73Zj5rEKmnq5YViGEynpKlQdWsje/s7SG1a5S7kaKWbSRmoUj3NlCZh3j33sPnTfEvf1FosXq0vhoLODDMqOyaDMXS1TwHApgASn3dvggsxz1UXmbsROtYXg/EiQqxVzxDnwrqIBvlOitbYw6KOpAWfXmr0bbo0WB0nKUizJQ+oCOA72olfFVAe8iXfRUdU+32UPv6HZrfMqKGtq3qpBmCJmlRe0abLZQqr+VZHUba+w21A3bLCDbi9WTDeqJZNFmnp7Bhz1iMYZqn0O3wT7AKblCJmMT8GyOtN6alTuHmxaEG2PmH6nxbWTy1pxuu+TBH6rwHYzuBYoZtkT2Vi4npeI62xb66O1dNquGL335RdzMoRsBqQMvfeJ9bo+m3KkPqH3Pe7c3mYeSSdFSOSvja6PnR//Al7jxh2b8X7/6Ke684WhvmP/c3lDCc8rkbcfWtyLtptDeEIqZUt9vcW1AumBi7YPV1w7ZX3WNzufE2X2bKrC9YYUQ6T151NzjJ9B5ZpUhO+qV116/7Eu5egy+7fCzev/wr/dW/D4wWNwQrEmaqlleTZlRnSJHPtsqtn8aPfgWfCNM25K9tkbalDsMiFt+ocxCyX5fcaOc89xoD4BpX9pWT9HSR89+X+NaoWhsHtAqEqxzhaR9WonDtIF4+HUs/NRjamlXs5+OrjxOIVvapyS3oOHkKPJJiFt2qigKG6mRJrlTFraEDgE36xi/PSPWBc3tGhRufLW1IdGbnn4kHDx/OKI7+bZS34d9ucN+pYwObCJBeC5SVIG33rvJm+E22jkQ5Xv+8Ft81+b7/L+/90fgOzXf98mv8mduf4l/+IUfIX51k9GOML9tXRq1gGpH8O9aJtShlxQxCxuHUZY26Y8+WD71QAhpdo+zfVsn1gTAhQ/d5x64NqLN3RevDkO0+CE8ywjy4slpjzIJl9VeTCHiWts2iZXgOrXb6ojeH+oPpSkBKXqLGxVTcI3ggt2nnSNUjjjkJPfmo876klmoiNOCes/x/myT97stmoOKya4ttUOdItTRsqB8k0ZlqjV8843lSBOBsKzqWXlDVl6uLlYaNjpk5X16xAkE10a0F5Fz/KwxZDydaAFO8GVPtLwaU3laRDSiTUqymEzgxuaicEBC2uMcBjWL9YBSJ8hmuSgWUIHRB2ophAWEWtj9Iyaw6kHagklbNKM3S2sQN7Im5zoJSBl5f3+DvabG73mKGbzz357nM1/6KC4K849EfGPL62IKxS6L+T9DY7fJez1b3wi2j9wF6xM1nS9SGOmDFQ4M7yktsW1Rse0t4NiB26dxbUR72TyL7VXllASAk1A1C6WDMFNmlIRhsHPy7VL+8SDcwXAvorXDNSgQoUgzeaxDY0pbiiYuJO2jio0CEaep2kiZzUuatsD1qTHcgeAbTxin2lwB9YoEa9NqJ00nj1gp3qxbLI0lpASRlV5Yx78RSaRRD+/RPgJZtGfEdfRtPyyaOfQ2Om28hRyKHqdui11HnDcmeu/R2RxtO9z2FnpnG+YdMp0hoYDC2XL3ICAKrovE3lkXxImw95LDt3D7i10ScEEsVzZ0UxJE85FIHAfcgUd6we96VDza1sQIxczSIMNICRVUu4KfmcjVQ31f2fh2YHbH096wJXwxD9b0baui2O+QabPch+37xVbPogdUWVj/45TxNEwKtD94uK72JLJoM0/Hh1mIxbweNZ+V5MvFgOJxhVW1LH6SX6tD5Dj9bTFPFitYVpQEh8RlozUXbI+02vfE0oZJL/KBfcoXHvzSnkUUtzgQpGfRrE0tPoVrzS8OI0uT9B1ISL5sL+lf2ztWJxw7ECCtJBava6idXa2Xdc7eo1MyoI6y1qIdgk9ADkBdZdKHV/uOYfCWlNWizYxZZAdDtUuM6NTGQ7q6wnc9urOLVBXxzk3ipMSNPIVTNt62JezsbmFLYU3Botbu70cCUdn+ikWnByHH0oRcTK287uDjSr8VGb3rKXehnKoJd2bPDzXsfaJg49uBG1+17hqLcSDDNY8rs6i9tYfVtrUtrY3JIiJOsHm2Ku7wfmt6Tx7F4q61aHPwaQ2Iat0ZFvuVZmGllGWN7DAGxNlSW8IwPT11fWg7S3ccj5HtTSTYRDyJSnlgPaFCKUjKSrJ623T6NBaz3NeF+NQv59m6LmVfDb5rD8VM0/1mwVFLl+zHFsH2H+xDlQJbJB9dBPUe8WEh5EWkXGzbacjFtvciml9/NC7wCPu1ay3azBXDWbeGRWvUGNCVIIsUxXLiuyo6b9B5g9sYW6VP1yN7M4sgj0fIeEwc1VAWeJ6ztqreWeXMLY8LsPHmHLxw8EJt6YtNND82tVKtdyLqYHbHE2rzb1HYejNQzCL7zxf0G8LkHeCdApeq5/qxVfAUs/RF0IOfW4vWuD227R2wfdmgiBcrXI9WD7vovhgCOrXVoI30XCZfLFYeQxF8vzKa/hSyaDNnhriU2SMKOpT4rPi8qW52KEsbloLalea/Bl12rXA235UiLWnH9dKCiaRZtEqx36ah0pWlMg4ZSikf2HXLgJSKCdYCW5Fqt8N9pLBl8ixlWg2r3SoFoZrk4wYW+YNaeqSPqcuGIiGtHkSXwtSVwoDVCPLq1k6yuLaB/egR5CzazJmxiCaLW4ozCVCcLCPJzi37IDnzOePOri2Th+dE6/Qg3TBt3acujbbFMrof8POI+2AXvGO0XdOPPN2mwwUYvzPH9ZH2RkUsHeN7FukNlVlgUQi1pzyISHTM7gqhho23ldGDwMFzlnHlOqwTZHSEKLhWbXunC0jTI02LzuY4Z9ld2rbmj4OtKga6jnhcj+fB4j4GWbSZs0VtqYhPIhss7urWz8oUdYoCuhadNckSH8kQSgOal36hWTM/jxZR7jroBTcPeBG6LY86Xe6dblcgtmUE4MZ+sc+rhc0C8o1FgcPIfNhiFnF9mjjfK65RnLc9YBdSLnRIFrTrYTZHvYcqQtPYPq33h5qzqXMnW9LH2KOFLNrMRRAD2kbLAnIpojosdUNc1NgOW0I4MavV9/D+PaSq0JtbaOoVhQjFNCBRCS/cXZzGtYFyz/qlzl7csDtTMn+100KMHGxvWnlfiEhvCRmuh/F7lrPsWwiVo96NlFNh9H5HtdMuBkq7PhLHqY9xjLA/tS4ckwlSjFdmFA2jLtNYkMGXFVtxHGpanpDUH/laR48za4Sa1dWIWeDVActUuNWlpDj7SY293UQhbloHi9qbr9rYBz5s2t9Jb3WsvonEUmhv2FSCcj/i+ojMOvM9hcUgLUmzaEWgOoiomO+qBRSzaJlWBz1uaoJ3XQqEeWezL4d95xAsUjx0YITl1o0IuOqwLysOcfFw5uJjVPqspWhzccDVRtN+5Go52mpm1NEPp/YdcRrNzx2Plh/8lDWkfY+88x5S11QHW1AWhM0aLRxaWEDKzXuzmq21ryn3zVoNOcFufwqqjL/dUI081YMG6QLt3Qn92COtbe241jos+nmPNAGtPd2tsQm/t64UQ6mgFg62NynK0jpT3LsP3i8mCTCbHVvsf6wP+xi+7VqKNu/PXnEemk3jDvm0xz1/0VrVrYo7dXfoe8J0aoO7NNp+7aSyNlEiCBbBtfm1dm53kI4xFCCkTojFfotrPW5ninQ9bnuEVA4XhkkGFmTyBy00LaHcJNY+ZUNZ1pb00QTrHDqyrpHSdYTpFDcawWhkBREnvean9G3XUrRXmWexcOBEVvy3odPgoYeL0j7kw+1RjVQl2jToflJdSsBwW1u2/AwRnTf4D/Zw3lF42xqKdQlOcTsH9iWQtogWk9lvbqHeE0cF6h1xc4zEiJ91Zo1DXAytJsQ0PCvi7x/gd2aLgnbpg3WkgBQo69EQ0aZZ+uXJJ5WiYGgXu4ioH+PLPi5ZtGfMdSwceGIGCxv7E63IoW2RolgUIMS5lbdJVS0GdaHRBNG16D6LbCIpS3j+lp2i621pXlfDAB8A4qhCa08s7Esgjgok2IQ715kQzYoG2yserq+ZWilh2tKhTx0VY+quuNomZrWdDizbyDJcq0c4YQzoY5BFmzk/NB5fKzr4uBqJs7lFj4fBW/v7aIi4NE1PNjcX+bo2NnJZRA6YmLwzCwlmYVcKzenNt3bTBlqPB3t88JvbztIm6xKtSkvuV0Wazr4AFq8lWdXZDJ03dr1lsfwyCiFZ0SN+vCra9YuKprOwuFm0mfNDV/ZpVxkscAho1wI2nEu7frnHWVXIaIRuTUxAs9lSEMNhht7CIhYgcqldzaHzp35TbWe1u8PM2DLN1hlujyq09JbZNDRl61YPZdla2nbEuSVTSF0jPrC4qph+X00kiTGlc6b7zsDiZtFmzo/Bpx0iyas+7qoFXqlsWd2v1fkc3jOxadsdPq44lAIZeRu9sWEW1t8/WKYNiqCbE0vOGAoQ2g7tU8G6iJ0jlQH6KvVaVrWEiaZF6soqkdrOunEMI0yGfOFwZDURA9qRvpCsVaobjZbPG3iKNXIW7TmRZ9iy4tMmi7vq464sDTWVrYk/XDSvbUec7ix828VhRUyYw3S6whNr81FpO7TrFhPm460ttHS4eWfWvA/mE6cWrjqbo31vacVdEm1UdD5H2xZfV2iZjt0HE/pqtlOypIeISysqRWHtdPo+rSqenrUS7Trtz+YZtjzs06bb4gTKeplckayTDB0Mh9pTkYVVc3W9aLm6QCT1S+4oBivmHUKZRNfh9me2JTM1i4p3qKst3TAtd1ErkZOisFlDhcftpfrbEJH9qWU0JXFKUVgyRdOkL4+TE/41BOJs/lSW9ShrJdq8P7tmHPVph9u+wtW1JdcPkdcYUC0Q5y3vcFjiDsvRul62XAWzpiHY8jYq7O2b4Le3bHk9U+iXTdb0wIrqZTxGioI4naXzd8vAlHPouEarYhlF7nu0OZJumJb8cVjSH8ruevg9OCsLO7BWos1cDzQEYtMsfbzVFL6mSe1oUhS3KJb1pmpJFEMkl9VoLZhw9vbtuW1nx2gaO3YKAGnaZyWJ3k0mdvy2Q9sO6XukrtCD2fJaUsuYY0dQRl3Z1vHLKPFjTAx4XLJoMxdPDIetlzizVnCofG3h4676uWmEpAnpsIg0RJgfHBZX0y4ttDjobOsopmWxG4+Q0Yj4/gcWFW7bRerkQoxYZPuQ7+p8atSW9mzTbRlqhldXEWdMFm3m8ln1dVPgahCMFMVin1a8h7pe/p2Twz2FNaJghQGD4EKw5IjimI+6qnWVaLvURNzb8nlUo7t7p0d7T/DXF83cTmqdegZk0WYun6ECSJdJF8QAaW6spmofKStctRyBiTikrg41VyMk4Q4tbYJFfM1vffjUcShYLwqzlJMxOq6Rg+npyQ8n+OsaWbaaOSeyaC+AnI98hKNlaCuR18HiSlktB3OlyhkAnc2t84WIjcxMPq42bRJ7Ok6Ih1MIZ6nr/1BG5wTFL4SpwfKJdW8fmub4LhPHvIahxc5imXwBZNFeADkf+QgrPiywDNoM1ssVuPFKIUFRwGgEXUvY2QVxuI2JLYdTHnDY2TUfdXMT4NCoEo0PF5YPllXbuMic0r4n7O4+1mtY9WHjPIs2s66cYkkXt1Ply0MtRFNp3qGqIOdMoCEurJt1hZBlm9V0jsXfiQOO2Tsdru1onetw/0nW8rgi9UW+8Qm9ik/Zv30asmgzZ85DEd+uf8j/075PWyQrH8E0OoMQrMpn9f4hMJWSL7RpLM939XmApmWtlBW2xXP6oOuj9y8s70NPOLw6GPKk7ffj3gQ5ff/2Kbiyol3NfroufuCz4tsO/uHyjmOygYY9zeMQd7gYHg5FYzXqyVZssIiLPMLDFlejIhwJID0KGg9HhB8lw2mwxGfMlRXtkP30ymuvH8rhXWeeGd92sKQnkfY0T35cDlvgo9ZKxIZY4R/KNrL7Wdz/kMWNT1hd82Gv6Wmf/xhcWdEOXGeLlDmGcDjpgqEB+iP6hxqCRZ+HfdnFPupjWrxFg7aHLfmhyqVL4NFHdWUy502yTsuf5Cw6WURqP5SYZt7WNTIeWxngE1i8oab36DJdqspSH1c7blwwV97SZq4hp0WPj1q41eeJO+SznugvDo3iRB65uuaor7u4ffTvh6X6OWY8fRhZtJmL57TocVUd2rddDGSuquUg5keIyC6iy4/KUV/3BN93WAVcJlm0l0AukE8ki6tHK3VWbi9806PR40dlWN5eULbSRXDlRLtOhe5PSi6QTxyNsB65LUWBjMe2b5tyhB8LEcusco74YbnEa8SVE20udH9GGZL6V9Cohyt2jj4fTo/gDrWwqmfaOeKyuRKiHawrrP9ebOYJcH5R3XPI8sYTLOywT3v0+ccQj2RMXQeuhGifZev6rGRJnYpGUHl0X/WaWc7HJe/TXjK/9uoPL76wXn71F3nltdcv+YougcGXfQyf8ypEcS+LK2FpM89AiuMxPuulHmeNuVTRPguR4gwLn/WpC8Ufw5e9zlyqaJ9lXzaTeVIuRbTZwp7MtUy8iMFGajwt59BDeB2RY3u5Ljkz5+Hots61+DCeM6fUFJ9QiPpk/Jj7S8+2k3hF+eX4C8f+P5+baFdFClmoT8PR9/Lrr306i/YZ4CTRnsny+OiHCkyk2V89G/KXXWaVD7O0mUzmipGTKzKZNSOLNpNZM7JoM5k1I4s2k1kzTo0ei8h/AO4+5TnuAu8/5THOgut0He+r6o+fxcVk1o9zjx6LyG+o6ifP9ST5OjLPEHl5nMmsGVm0mcyacRGi/acXcI5HIV9H5lqQM6IymTUjL48zmTUjizaTWTPOVbQi8ldF5H+IyBdE5HMi8n3neb4TruHHReRLIvIVEXn1os+fruElEfkVEfmiiPyOiPyty7iOzPXgXH1aEflTwO+q6n0R+XPAz6rqnzy3Ez58fg98Gfgx4E3g88BPq+oXL+oa0nW8ALygqm+IyBbwm8BPXfR1ZK4H52ppVfVzqno/3fyvwMfP83zH8EPAV1T191W1Bf4l8JMXfA2o6juq+kb6fQ/4XeDFi76OzPXgIn3avw780gWeD0wY31q5/SaXLBYReRn4AeDXL/M6MuvLhTR2E5E/i4n2UxdxvquKiGwC/xr4GVXdvezryawnZ25pReRvishvpZ+PicgfBz4D/KSqfnDW5/sQ3gJeWrn98XTfhSMiJSbYf6Gq/+YyriFzPTjvQNQngNeBv6aqnzu3E518/gILRP0IJtbPA39FVX/ngq9DgH8O3FPVn7nIc2euH+ct2s8AfwH4Rrqrv+gKFxH5CeAfAx74eVX9Bxd5/nQNnwL+M/AFYJgc9XdU9d9f9LWcQE6Lu3qc2HEzpzFmIIv2KnKiaHNGVCazZmTRZjJrRhZtJrNm5Pm0mWeCU+YirR3Z0maeCVbHqr786i/yymuvX/IVPTlrL1oR+VkR+dvp978vIj/6FMf6eRF5T0R+++yuMHOV+LVXf5ivv/bph2ZPrRNrL9pVVPXvqep/eopDfBbIrUkzV5q1FK2I/F0R+bKI/Bfgj67c/1kR+Yvp96+LyM+ldMrfEJEfFJH/KCJfFZG/cdxxVfVXgXsX8yoymSdj7QJRIvIngL8MfD92/W9g9anH8U1V/X4R+UeYFX0FGAG/DfyTc7/YzFqwbrOU1060wJ8G/q2qTgFE5N+d8tzhsS8Am6mWdU9EGhG5qaoPzvdSM+vAapAKLFB1lVlH0T4OTfo3rvw+3L7urz1zCi/eHC/E+eLN8SVfzeOxjh/cXwU+KyI/h13/nwf+z8u9pMxVZXV/dpWrvPz9MNZOtKnP0r8C/jvwHlZudyaIyP8N/Bngroi8CfwfqvrPzur4mYvn6NL3OrB2ogVI5XUPldip6v+28vvLK79/FgtEPfTYkb//6TO7yEzmnFjLLZ9M5llmLS1tJnMaq1s46xZkehSyaDPXjuvox66Sl8eZzJqRRZvJrBl5eZzJHOFo4sVV29PNos1kjrAq0quY0piXx5nMmpFFm8msGVm0mcyakX3azLXhpOKA60YWbebacN2TKgby8jiTWTOyaDOZNSOLNpNZM7JoM5k1I4s2k1kzsmgzmTUjizaTWTPyPm1m7TnPpIqh4ucqVftk0WbWnvNMqhiEepWqffLyOJNZM7JoM5k1I4s2k1kzsmgzmTUjizaTWTOyaDOZNSOLNpNZM7JoM5lHYEiyeOW11y/7UnJyRWY9ueh5PVcpySKLNrOWPCutZY4jL48zmTUjizaTWTOyaDOZNSOLNpNZM3IgKrNWPCsNyU8jizazVjzLUeOBvDzOZNaMLNpM5jG4CplReXmcyTwGVyEzKlvaTGbNyKLNZNaMLNpMZs3Ios1k1ows2kxmzciizWTWjLzlk8k8AcN+7fD7RY4MyaLNZJ6AVZFe9J5tXh5nMmtGFm0ms2bk5XFmLcgleUuyaDNrQS7JW5KXx5nMmpFFm8msGVm0mcyakX3azJUmB6AeJos2c6XJAaiHycvjTGbNyKLNZJ6Si+4blZfHmcxTctF9o7KlzWTWjCzaTGbNyKLNZNaMLNpM5oy4qIBUDkRlMmfERQWksqXNZNaMLNpMZs3Ios1k1ozs02auHEORAJALBY4hizZz5chFAqeTRZu5MuQyvEcjizZzZbguFva8G5ln0WYyZ8x5NzLP0eNM5hw5jyypbGkzl8519mXPI0sqizZz6VwXX/aiyMvjTGbNyKLNZNaMvDzOXAo56+nJyaLNXArZj31y8vI4k7kAznLrJ1vazIVynbd3TuMst36yaDMXSl4WPz1ZtJkL4Vm1sEdZzUteve9x8pNFVc/6ujLrx7l8CI5GiM86cf66sPqFtvIeyUnPz6LNwFOK9uiH7oQPYeZDWP2S+/prn86izWSuC3nLJ5NZM7JoM5k1I4s2k1kzsmgzmTUj79NmEJH/ANx9ysPcBd4/g8t5Wq7CdZzFNbyvqj9+3AM5epw5E0TkN1T1k/k6zv8a8vI4k1kzsmgzmTUjizZzVvzTy76AxFW4jnO9huzTZjJrRra0mcyakUWbeSpE5K+KyP8QkS+IyOdE5Psu6Tp+XES+JCJfEZFXL+kaXhKRXxGRL4rI74jI3zqX8+TlceZpEJE/Bfyuqt4XkT8H/Kyq/skLvgYPfBn4MeBN4PPAT6vqFy/4Ol4AXlDVN0RkC/hN4KfO+jqypc08Far6OVW9n27+V+Djl3AZPwR8RVV/X1Vb4F8CP3nRF6Gq76jqG+n3PeB3gRfP+jxZtJmz5K8Dv3QJ530R+NbK7Tc5B7E8DiLyMvADwK+f9bFzGmPmTBCRP4uJ9lOXfS2XjYhsAv8a+BlV3T3r42dLm3lsRORvishvpZ+PicgfBz4D/KSqfnAJl/QW8NLK7Y+n+y4cESkxwf4LVf0353KOHIjKPA0i8gngdeCvqernLukaCiwQ9SOYWD8P/BVV/Z0Lvg4B/jlwT1V/5tzOk0WbeRpE5DPAXwC+ke7qLyNhX0R+AvjHgAd+XlX/wSVcw6eA/wx8AYjp7r+jqv/+TM+TRZvJrBfZp81k1ows2kxmzciizWTWjCzaTGbNyKLNZNaMLNrMlUdEflZE/nb6/e+LyI8+4XEupArnvMlpjJm1QlX/3lP8eQ/876tVOCLyyxddDfS0ZEubuZKIyN8VkS+LyH8B/ujK/Z8Vkb+Yfv+6iPxcSqf8DRH5QRH5jyLyVRH5G0ePeVFVOOdNtrSZK4eI/AngLwPfj31G38BqU4/jm6r6/SLyj4DPAq8AI+C3gX9yyjle5pyqcM6bLNrMVeRPA/9WVacAIvLvTnnu8NgXgM1kQfdEpBGRm6r64OgfnHcVznmTl8eZdadJ/8aV34fbDxmli6jCOW+yaDNXkV8FfkpExilg9OfP4qCpCuefYe1x/uFZHPMyyKLNXDlSsOhfAf8d64Tx+TM69CvA/wr88Eo98E+c0bEvjFzlk8msGdnSZjJrRhZtJrNmZNFmMmtGFm0ms2Zk0WYya0YWbSazZmTRZjJrRhZtJrNm/P9CLm4wL/KaXAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "_ = pairplot(posterior.sample((10_000,)), limits=[[-3, 3], [-3, 3]], figsize=(4, 4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "sensitivity = ActiveSubspace(posterior)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This time, we begin by drawing samples from the posterior and then computing the desired property for each of the samples (i.e. you will probably have to run simulations for each theta and extract the property from the simulation output). As an example, we assume that the property is just the cube of the first dimension of the simulation output:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "56a429406ff641d685f63dbce716fbe1",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Drawing 5000 posterior samples:   0%|          | 0/5000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d1b06137f0dd419693a31f4117134f96",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Running 5000 simulations.:   0%|          | 0/5000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "theta, x = simulate_for_sbi(simulator, posterior, 5000)\n",
    "property_ = x[:, :1] ** 3  # E.g. metabolic cost."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To investigate the sensitivity of a given parameter, we train a neural network to predict the `property_` from the parameters and then analyse the neural network as above:\n",
    "\n",
    "$M = \\mathbb{E}_{p(\\theta|x_o)}[\\nabla_{\\theta}f(\\theta)^T \\nabla_{\\theta}f(\\theta)$]  \n",
    "\n",
    "where $f(\\cdot)$ is the trained neural network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Training neural network. Epochs trained:  24"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "dee98f26b7c24c479efd0fdc888c9812",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Drawing 1000 posterior samples:   0%|          | 0/1000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "_ = sensitivity.add_property(theta, property_).train()\n",
    "e_vals, e_vecs = sensitivity.find_directions()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenvalues: \n",
      " tensor([2.8801e-06, 6.1131e-05]) \n",
      "\n",
      "Eigenvectors: \n",
      " tensor([[ 0.0362,  0.9993],\n",
      "        [ 0.9993, -0.0362]])\n"
     ]
    }
   ],
   "source": [
    "print(\"Eigenvalues: \\n\", e_vals, \"\\n\")\n",
    "print(\"Eigenvectors: \\n\", e_vecs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, one of the eigenvalues is much smaller than the other one. The larger eigenvalue represents (approximately) the vector `[1.0, 0.0]`. This makes sense, because only the `property_` is influenced only by the first output which, in turn, is influenced only by the first parameter."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
